rm(list = ls())

library(dplyr)
library(magrittr)
library(here)
library(fixest)

# Define relative path to save in the script's directory
file_path <- here()

# You should download "expcor_cls_apsr.rda" to run this script: Tai, Yuehong 'Cassandra'; Hu, Yue; Solt, Frederick, 2022, "Replication Data for: Democracy, Public Support, and Measurement Uncertainty", https://doi.org/10.7910/DVN/XAUF3H 


load(here('expcor_cls_apsr.rda'))


# Adding V-DEM vars  ----

## Adding vdem variables to further breakdown the liberal component (note that you need "vdemdata" library to run this:)

vdem <- vdemdata::vdem

#e_migdppcln missing
vdemvars <- vdem %>% select(country_name, year, v2x_polyarchy, v2x_liberal, v2xcl_rol, v2x_jucon, v2xlg_legcon, v2clrspct, v2xnp_regcorr,v2x_execorr, v2x_corr, v2xpe_exlecon, v2xpe_exlgender, v2xpe_exlgeo, v2xpe_exlpol, v2xpe_exlsocgr, v2peapsecon, v2peapsgen, v2peapspol, v2peapssoc, v2peasjsoecon, v2pepwrses, v2clacjust, v2x_libdem, v2x_polyarchy, v2jureform, v2jureform_ord, v2jupoatck, v2jureview, e_gdp, e_wbgi_gee, v2exbribe, v2exembez, v2x_pubcorr, e_wbgi_cce, e_ti_cpi, v2x_horacc, v2x_veracc, v2x_diagacc, e_miinflat, v2x_accountability, e_total_fuel_income_pc, e_miurbpop, e_civil_war, e_pt_coup, e_peaveduc, v2x_clphy, v2caviol, e_total_oil_income_pc)

colnames(vdemvars)[which(colnames(vdemvars)=='country_name')] <- 'Country'
colnames(vdemvars)[which(colnames(vdemvars)=='year')] <- 'Year'

## Fixing some country names before merging

vdemvars$Country <- ifelse(vdemvars$Country=='Trinidad and Tobago', 'Trinidad & Tobago', vdemvars$Country)

vdemvars$Country <- ifelse(vdemvars$Country=='United States of America', 'United States', vdemvars$Country)

vdemvars$Country <- ifelse(vdemvars$Country=='Burma/Myanmar', 'Myanmar', vdemvars$Country)

vdemvars$Country <- ifelse(vdemvars$Country=='Palestine/West Bank', 'Palestine', vdemvars$Country)

vdemvars$Country <- ifelse(vdemvars$Country=='Bosnia and Herzegovina', 'Bosnia & Herzegovina', vdemvars$Country)

vdemvars$Country <- ifelse(vdemvars$Country=='Ivory Coast', 'Côte d’Ivoire', vdemvars$Country)

vdemvars$Country <- ifelse(vdemvars$Country=='Türkiye', 'Turkey', vdemvars$Country)


## Add Swiid data (Version 9.7). It is available here: https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/LM4OWF


swiid <- read.csv(here('swiid9_7_summary.csv'))

swiid <- swiid %>% select(country, year, gini_disp, gini_disp_se)

# fixing some country names before merging
swiid$country <- ifelse(swiid$country=="Korea", 'South Korea', swiid$country)


swiid$country <- ifelse(swiid$country=='Palestinian Territories', 'Palestine', swiid$country)

swiid$country <- ifelse(swiid$country=='São Tomé and Príncipe', 'Sao Tome and Principe', swiid$country)

swiid$country <- ifelse(swiid$country=='Bosnia and Herzegovina', 'Bosnia & Herzegovina', swiid$country)

swiid$country <- ifelse(swiid$country=="Côte d'Ivoire", 'Côte d’Ivoire', swiid$country)

swiid$country <- ifelse(swiid$country=="Czech Republic", 'Czechia', swiid$country)

swiid$country <- ifelse(swiid$country=="Trinidad and Tobago", 'Trinidad & Tobago', swiid$country)


swiid <- swiid %>% rename(Country = country, Year = year)

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

# Independent variable: Judicial constraints ----

# Define a function to process each dataset
process_dataset <- function(data, vdemvars, swiid) {
  
  # Merge data
  data <- data %>% 
    rename(Country = country, Year = year) %>%
    left_join(vdemvars, by = c("Country", "Year")) %>%
    left_join(swiid, by = c("Country", "Year")) %>%
    group_by(Country) %>%
    mutate(
      e_wbgi_cce_imp = zoo::na.approx(e_wbgi_cce, na.rm = FALSE),
      gdp_change_pp = ((e_gdp * 100) / lag(e_gdp)) - 100
    ) %>%
    ungroup()
  
  # Define the models
  models <- list(
    feols(SupDem_trim ~ v2x_jucon | Country, data = data, vcov = ~Country),
    feols(SupDem_trim ~ v2x_jucon | Country + Year, data = data, vcov = ~Country),
    feols(SupDem_trim ~ v2x_jucon + gdp_change_pp + gini_disp | Country, data = data, vcov = ~Country),
    feols(SupDem_trim ~ v2x_jucon + gdp_change_pp + gini_disp | Country + Year, data = data, vcov = ~Country),
    feols(SupDem_trim ~ v2x_jucon + gdp_change_pp + gini_disp + e_peaveduc | Country, data = data, vcov = ~Country),
    feols(SupDem_trim ~ v2x_jucon + gdp_change_pp + gini_disp + e_peaveduc | Country + Year, data = data, vcov = ~Country),
    feols(SupDem_trim ~ v2x_jucon + gdp_change_pp + gini_disp + e_peaveduc + e_wbgi_cce_imp | Country, data = data, vcov = ~Country),
    feols(SupDem_trim ~ v2x_jucon + gdp_change_pp + gini_disp + e_peaveduc + e_wbgi_cce_imp | Country + Year, data = data, vcov = ~Country),
    feols(SupDem_trim ~ v2x_jucon + gdp_change_pp + gini_disp + e_peaveduc + e_wbgi_cce_imp + v2x_clphy | Country, data = data, vcov = ~Country),
    feols(SupDem_trim ~ v2x_jucon + gdp_change_pp + gini_disp + e_peaveduc + e_wbgi_cce_imp + v2x_clphy | Country + Year, data = data, vcov = ~Country)
  )
  
  # Extract model summaries
  model_summaries <- lapply(models, function(model) {
    # Extract coefficients and model info
    tidy_model <- broom::tidy(model)
    tidy_model$model <- deparse(substitute(model))
    tidy_model
  })
  
  # Combine model summaries into one data frame
  bind_rows(model_summaries, .id = "model_id")
}

# Process all datasets in the list
all_results <- lapply(expcor_cls_apsr, function(dataset) {
  tryCatch(
    process_dataset(dataset, vdemvars, swiid),
    error = function(e) {
      message("Error processing dataset: ", e)
      NULL
    }
  )
})

# Combine all results into a single data frame
final_results <- bind_rows(all_results, .id = "dataset_id")


# Save results for later use
saveRDS(final_results, file = here('regression_results_judicial_const_uncert.rds'))


# Legislative constraints ----

# Define a function to process each dataset
process_dataset_leg <- function(data, vdemvars, swiid) {
  
  # Merge data
  data <- data %>% 
    rename(Country = country, Year = year) %>%
    left_join(vdemvars, by = c("Country", "Year")) %>%
    left_join(swiid, by = c("Country", "Year")) %>%
    group_by(Country) %>%
    mutate(
      e_wbgi_cce_imp = zoo::na.approx(e_wbgi_cce, na.rm = FALSE),
      gdp_change_pp = ((e_gdp * 100) / lag(e_gdp)) - 100
    ) %>%
    ungroup()
  
  # Define the models
  models <- list(
    feols(SupDem_trim ~ v2xlg_legcon | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xlg_legcon | Country + Year, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xlg_legcon + gdp_change_pp + gini_disp | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xlg_legcon + gdp_change_pp + gini_disp | Country + Year, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xlg_legcon + gdp_change_pp + gini_disp + e_peaveduc | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xlg_legcon + gdp_change_pp + gini_disp + e_peaveduc | Country + Year, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xlg_legcon + gdp_change_pp + gini_disp + e_peaveduc + e_wbgi_cce_imp | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xlg_legcon + gdp_change_pp + gini_disp + e_peaveduc + e_wbgi_cce_imp | Country + Year, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xlg_legcon + gdp_change_pp + gini_disp + e_peaveduc + e_wbgi_cce_imp + v2x_clphy | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xlg_legcon + gdp_change_pp + gini_disp + e_peaveduc + e_wbgi_cce_imp + v2x_clphy | Country + Year, data = data, vcov = ~Country)
    
  )
  
  # Extract model summaries
  model_summaries <- lapply(models, function(model) {
    # Extract coefficients and model info
    tidy_model <- broom::tidy(model)
    tidy_model$model <- deparse(substitute(model))
    tidy_model
  })
  
  # Combine model summaries into one data frame
  bind_rows(model_summaries, .id = "model_id")
}

# Process all datasets in the list
all_results_leg <- lapply(expcor_cls_apsr, function(dataset) {
  tryCatch(
    process_dataset_leg(dataset, vdemvars, swiid),
    error = function(e) {
      message("Error processing dataset: ", e)
      NULL
    }
  )
})

# Combine all results into a single data frame
final_results_leg <- bind_rows(all_results_leg, .id = "dataset_id")

# Save results for later use
saveRDS(final_results_leg, file = here('regression_results_legislative_const_uncert.rds'))



# Ind liberties ----

# Define a function to process each dataset
process_dataset_indlib <- function(data, vdemvars, swiid) {
  
  # Merge data
  data <- data %>% 
    rename(Country = country, Year = year) %>%
    left_join(vdemvars, by = c("Country", "Year")) %>%
    left_join(swiid, by = c("Country", "Year")) %>%
    group_by(Country) %>%
    mutate(
      e_wbgi_cce_imp = zoo::na.approx(e_wbgi_cce, na.rm = FALSE),
      gdp_change_pp = ((e_gdp * 100) / lag(e_gdp)) - 100
    ) %>%
    ungroup()
  
  # Define the models
  models <- list(
    feols(SupDem_trim ~ v2xcl_rol | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xcl_rol | Country + Year, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xcl_rol + gdp_change_pp + gini_disp | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xcl_rol + gdp_change_pp + gini_disp | Country + Year, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xcl_rol + gdp_change_pp + gini_disp + e_peaveduc | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xcl_rol + gdp_change_pp + gini_disp + e_peaveduc | Country + Year, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xcl_rol + gdp_change_pp + gini_disp + e_peaveduc + e_wbgi_cce_imp | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xcl_rol + gdp_change_pp + gini_disp + e_peaveduc + e_wbgi_cce_imp | Country + Year, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xcl_rol + gdp_change_pp + gini_disp + e_peaveduc + e_wbgi_cce_imp + v2x_clphy | Country, data = data, vcov = ~Country),
    
    feols(SupDem_trim ~ v2xcl_rol + gdp_change_pp + gini_disp + e_peaveduc + e_wbgi_cce_imp + v2x_clphy | Country + Year, data = data, vcov = ~Country)
    
  )
  
  # Extract model summaries
  model_summaries <- lapply(models, function(model) {
    # Extract coefficients and model info
    tidy_model <- broom::tidy(model)
    tidy_model$model <- deparse(substitute(model))
    tidy_model
  })
  
  # Combine model summaries into one data frame
  bind_rows(model_summaries, .id = "model_id")
}

# Process all datasets in the list
all_results_indlib <- lapply(expcor_cls_apsr, function(dataset) {
  tryCatch(
    process_dataset_indlib(dataset, vdemvars, swiid),
    error = function(e) {
      message("Error processing dataset: ", e)
      NULL
    }
  )
})

# Combine all results into a single data frame
final_results_indlib <- bind_rows(all_results_indlib, .id = "dataset_id")

# Save results for later use
saveRDS(final_results_indlib, file = here('regression_results_indliberties_const_uncert.rds'))


